!> 读取节点和单元数据的模块
module mesh_reader_mod
    use iso_fortran_env, only: real64
    implicit none
contains
    ! 读取矩阵数据的函数
    ! 输入参数:
    !   filename - 文件名(字符串)
    ! 输出参数:
    !   matrix   - 动态分配的实数矩阵
    !   nrows    - 矩阵行数
    !   ncols    - 矩阵列数
    !   ierr     - 错误代码(0=成功, 1=文件打开失败, 2=读取失败)
    ! 文件格式:
    !    # 这是一个矩阵数据
    !    # rows=3, cols=3
    !    1.0  2.0  3.0
    !    4.0  5.0  6.0
    !    7.0  8.0  9.0
    subroutine loadtxt(filename, matrix, nrows, ncols, ierr)
        character(len=*), intent(in) :: filename
        real, allocatable, intent(out) :: matrix(:,:)
        integer, intent(out) :: nrows, ncols, ierr
        
        integer :: fileunit, i, stat
        character(len=256) :: line
        logical :: size_found
        
        ! 初始化错误代码
        ierr = 0
        
        ! 打开文件
        open(newunit=fileunit, file=filename, status='old', action='read', iostat=stat)
        if (stat /= 0) then
            ierr = 1  ! 文件打开失败
            return
        end if
        
        ! 初始化行列数
        nrows = 0
        ncols = 0
        size_found = .false.
        
        ! 读取文件头信息
        read(fileunit, *, iostat=stat) nrows,ncols
        if (stat /= 0) then
            ierr = 0  ! 不读取行列数
            return
        end if 
        ! 如果没有找到行列数信息，则通过数据统计确定
        if (.not. size_found) then
            ! 先统计行数
            nrows = 0
            do
                read(fileunit, *, iostat=stat)
                if (stat /= 0) exit
                nrows = nrows + 1
            end do
            
            ! 统计列数(回到文件开头)
            rewind(fileunit)
            do
                read(fileunit, '(A)', iostat=stat) line
                if (stat /= 0) exit
                if (len_trim(line) == 0 .or. line(1:1) == '#') cycle
                exit
            end do
            backspace(fileunit)
            
            ! 读取第一行数据确定列数
            ncols = 0
            do
                read(fileunit, *, iostat=stat)
                if (stat /= 0) exit
                ncols = ncols + 1
                if (ncols > 100) exit  ! 防止无限循环
            end do
            
            ! 重置文件指针
            rewind(fileunit)
            ! 跳过注释行
            do
                read(fileunit, '(A)', iostat=stat) line
                if (stat /= 0) exit
                if (len_trim(line) == 0 .or. line(1:1) == '#') cycle
                exit
            end do
            backspace(fileunit)
        end if
        
        ! 检查行列数是否有效
        if (nrows <= 0 .or. ncols <= 0) then
            ierr = 2  ! 无效的行列数
            close(fileunit)
            return
        end if
        
        ! 分配矩阵内存
        allocate(matrix(nrows, ncols), stat=stat)
        if (stat /= 0) then
            ierr = 3  ! 内存分配失败
            close(fileunit)
            return
        end if
        
        ! 读取矩阵数据
        do i = 1, nrows
            read(fileunit, *, iostat=stat) matrix(i, :)
            if (stat /= 0) then
                ierr = 3  ! 数据读取失败
                deallocate(matrix)
                exit
            end if
        end do
        
        ! 关闭文件
        close(fileunit)
    end subroutine loadtxt

    !> 加载节点坐标数据,coords(节点数,3):节点坐标数组; nnum:节点数
    subroutine load_nodes(filename,coords,nnum,ierr)
        character(len=*), intent(in) :: filename
        real(real64), allocatable, intent(out) :: coords(:,:)
        integer, intent(out) :: nnum,ierr
        integer :: fileunit,i,stat,nlabel
        character(len=200) :: iomsg
        
        ierr = 0
        open(newunit=fileunit, file=filename, status='old', action='read', iostat=stat,iomsg=iomsg)
        if (stat /= 0) then
            ierr = 1  ! 文件打开失败
            print*,"Error in : file open failed! ",iomsg
            return
        end if
        
        ! 读取节点数
        read(fileunit,*,iostat=stat) nnum
        
        ! 分配坐标数组内存
        allocate(coords(nnum,3),stat=stat)
        if (stat /= 0) then
            ierr = 3  ! 内存分配失败
            close(fileunit)
            return
        end if
        
        ! 逐行读取坐标数据
        do i=1,nnum
            read(fileunit,*,iostat=stat) nlabel,coords(i,1),coords(i,2),coords(i,3)
            if (stat /= 0) then
                ierr = 4  ! 数据读取失败
                deallocate(coords)
                close(fileunit)
                return
            end if 
        end do

        ! 关闭文件
        close(fileunit)
    end subroutine load_nodes
    
    !> 加载单元连接数据,elnodes(单元数,elnn):单元节点数组,elconn:单元节点数
    subroutine load_elements(filename,elnn,elnodes,ierr)
        character(len=*), intent(in) :: filename
        integer,intent(in)::elnn
        integer, allocatable, intent(out) :: elnodes(:,:)
        integer, intent(out) :: ierr
        integer :: fileunit,i,stat,elnum,elabel
        
        ierr = 0
        open(newunit=fileunit, file=filename, status='old',&
            action='read', iostat=stat)
        if (stat /= 0) then
            ierr = 1  ! 文件打开失败
            return
        end if
        
        ! 读取单元数
        read(fileunit,*,iostat=stat) elnum
        
        ! 分配单元节点数组内存
        allocate(elnodes(elnum,elnn),stat=stat)
        if (stat /= 0) then
            ierr = 2  ! 内存分配失败
            close(fileunit)
            return
        end if
        
        ! 逐行读取坐标数据
        do i=1,elnum
            read(fileunit,*,iostat=stat) elabel, elnodes(i,:)
            if (stat /= 0) then
                deallocate(elnodes)
                close(fileunit)
                ierr = 3  ! 数据读取失败
                write(*,'(A,I0,A)') "Error in : element read ",i," failed!"
                return
            end if 
        end do
        
        ! 关闭文件
        close(fileunit)
    end subroutine load_elements

end module mesh_reader_mod